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A compressed sensing scheme for near-field imaging of corrugations of relative 
sparse Fourier components is proposed. The scheme employs random sparse 
measurement of near field to recover the angular spectrum of the scattered 
field. It is shown heuristically and numerically that under the Rayleigh hy- 
pothesis the angular spectrum is compressible and amenable to compressed 
sensing techniques. 

Iteration schemes are developed for recovering the surface profile from the 
angular spectrum. The proposed nonlinear least squares in the Fourier basis 
produces accurate reconstructions even when the Rayleigh hypothesis is known 
to be false. 
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1. Introduction 

Rough surface scattering is of fundamental interest in optics, radiowave propagation and 
acoustics jHEl|28] and forms the basis of near-field imaging which is the operation principle 
behind such instruments as scanning near-field optical microscopy [31 [161 [201 [25] an d near 
field acoustic microscopy [TJJ]. Near- field imaging is a microscopic technique that breaks the 
diffraction limit by exploiting the properties of evanescent waves. The signal is collected by 
placing the detector in a distance much smaller than wavelength A to the specimen surface. 
An image of the surface is obtained by mechanically moving the probe in a raster scan of the 
specimen, line by line, and recording the probe-surface interaction as a function of position. 
This leads to long scan times for large sample areas or high resolution imaging. 

Typically near-field imaging is analyzed by assuming a continuum or dense set of data 
points [T5|l24p27j . In the present work, we focus on the setting of sparse, discrete measurement 
of near-field from the perspective of compressed sensing theory. This is an extension of the 
work [13] on potential scattering to the case of rough surface scattering. Surface scattering 
involves the geometry (i.e. topography) of scatterers and is technically more challenging to 
deal with than potential scattering. 



Consider the scattering problem for a corrugation profile described by the function z = 
h(x). For simplicity of presentation, we will focus on the case of two-dimensional scalar wave 
with the Dirichlet boundary condition. The total field u tot satisfies 

Au tot + k 2 u tot = in tt C M 2 , k > (1) 

u tot = on dn, (2) 

where 

fi = { r = (x,z) £ M 2 : z > h(x)}, h £ C(E) fl L°°(K.). (3) 

The total field models the sound pressure wave or electromagnetic waves in the TE-mode. 
The Dirichlet boundary condition corresponds to the sound-soft boundary condition in acous- 
tics and in electromagnetism the perfectly conducting boundary condition. Our approach can 
be easily extended to the three dimensional case as well as to the Neumann boundary con- 
dition, corresponding to acoustically hard obstacles, and the Robin boundary condition. 

As usual in scattering problem, we write u tot = u mc + u where both the scattered wave u 
and the incident wave u mc satisfy the Helmholtz equation. The Dirichlet condition becomes 
u = —u mc on dtt. 

In this paper, we focus on the case of periodic surfaces which include diffraction gratings, 
an important class ofoptical elements. We assume that h has period L and u mc is the plane 
incident wave 

u inc (r) = e lk2 r = e^cose-zsme) ? ^ = ( cos ^ ; _ sin ^) ; < 6 < tt. (4) 
Observe that on the boundary z = h{x) 

u(x + L, h(x)) = - e ik ( (x+L) cos e - h{x) sine ) = e ikLcos9 u{x, h(x)). (5) 
Hence we look for the (L, /c cos ^-quasi-periodic (or Floquet periodic) solution satisfying 

u{x + nL, z) = e [nLkcose u{x, z) for all (x, z) e tt, n £ Z. (6) 
In particular, if 9 = ir/2, then u is L-periodic. To fix the idea, we set L = 2ir. 

2. Radiation condition and Rayleigh hypothesis 

The existence and uniqueness can be proved under the quasi-periodicity and the radiation 
conditions on the solution u [18]. The well-posedness for general nonperiodic rough surfaces 
is given in [1]. Below we discuss the Fourier representation of the scattered field and the 
associated Rayleigh hypothesis. 

For x £ [— 7r, 7r), z > sup h = /i max , we write the scattered field as the Fourier series 

u(x,z) = Y^u n (z)e l{n+kcose)x = J2 u n(z)e ikanX (7) 

with 

Tl 

a n = - + cos6 (8) 
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Fig. 1. Surface topography. 



where u n satisfies 

u n + k 2 (l-a 2 n )u n = 0. (9) 

Solving Eq.(|nj) and imposing the boundedness of u n at z = oo we obtain the general solution 

as 

u(x, z) = J2\a n \<i a n e lk ( a " x ~P nZ } (incoming waves) 
+ Ekki ^qM^x+^z) (outgoing waves) 
+ ^2\ a i >:L c^e 1 ^ ™ 2 ^™^ (evanescent waves) (10) 



where /3 n is given by 




\ an \ z ! ^ 



The Rayleigh radiation condition for the region above the grooves z > h max amounts to 
dropping the incoming waves in Eq.( jT0i) : 

u(x, z) = JZ u n e [k ^ x+ ^ , z > /i max . (12) 

However, in the region inside the grooves z < h max multiple scattering may occur and Eq. (|12|) 
may not represent the true scattered wave in this region. For shallow corrugations, Eq. (|12|) 
should hold in the grooves and this is the Rayleigh hypothesis. For instance, the Rayleigh 
hypothesis holds for the sinusoidal profile h(x) = bsm(ax) with \ab\ < 0.448 [22l|23l[28]. On 
the other hand, for a general periodic surface the validity of the Rayleigh hypothesis may 
be difficult to assess [17]. The failure of the Rayleigh hypothesis Eq. (fl2]) manifests in the 
breaking down of the analytic continuation of Eq. (fT2l) inside grooves. 
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3. Inverse scattering formulation 

Inverse scattering seeks to reconstruct h(x) by transmitting incident wave u mc and measuring 
the scattered field u at certain locations. Moreover, in order to resolve subwavelength struc- 
ture which is hidden in the evanescent waves the measurement should be carried out in the 
near-field. 

Due to the quasi-periodicity, we may consider the scattered field u in the union of 

vl v = {(x,z)en, is [-7T, tt)} (13) 

and 

r = {( X ,Z) G X G [— TT, 7T)}. (14) 

For z > h(x) we have the outgoing scattered wave representation Eq. (jl2p with [2] 



ik{ anX > +M *>)) (_dv^p_ \ J l+ - h 2 {x , )dx ,, (15) 
\ ou' r'er / v 



du tot (r') 

u n = . , „ / e 



To ensure (3 n ^ in Eq.f fT5|) . we assume that \a n \ ^ 1, i.e., 

7^1, neZ (16) 



cos 6 H — 
k 



to avoid all grazing modes. In the case of normal incidence 9 = ir/2, Eq. (fT6]) means that the 
wavenumber k is not an integer. 

A key assumption for our approach is that h has a small number of significant Fourier 
coefficients, namely the Fourier coefficients are sparse or compressible. Writing h(x) = 
^2 neZ h n e mx , we say that h = {h n } is s-sparse if ||h||o, the number of nonzero elements 
of h, is less or equal to a small integer s. Note that h*_ n = h n since h is real-valued. Without 
loss of generality, we assume that h = 0. 

For reconstruction of h(x) we utilize the sparsity of h which surprisingly yields com- 
pressibility of the scattering amplitude u n . Compressive sensing techniques can then used to 
effectively recover those modes. 

Let (xj, zq), j = 1, 2, . . . , m, be the sensor locations for measuring the scattered field where 
zq > /i max is fixed and Xj are randomly and independently chosen from [—tt, it) according to 
the uniform distribution. 

In view of the identify 

u(xj, z )e- ikcos9x i = e [nxs u n e [k ^ zo (17) 

from Eq.( fT2l . let us consider the following inverse problem Y = AX, with entries 

X n = u n e ik ^ Z0 V^ (18) 
Yj = u(x j ,z )e- ikx ^ os9 (19) 

A = [A^ n } = ^=^ (20) 



where n is restricted to a finite, but sufficiently large interval ranged from —N/2 to N/2 — 1. 
In general the system Eq. fll8|) -Eq. fl20p is highly underdetermined for any m < oo. 

Surprisingly, sparse Fourier coefficients {h n } give rise to sparse or compressible {«„} and 
therefore X which can be reconstructed by compressed sensing. 
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4. Compressive sensing (CS) 

The main thrust of compressed sensing [71fT5] is to convert the noisy underdetermined system 

Y = AX + E (21) 

into the L 1 -based optimization problem 

minHXl^ subject to \\Y - AX\\ 2 < e = \\E\\ 2 (22) 

where E is the external noise vector. Eq.f l22p is called the Basis Pursuit (BP) [SJ. In addition 
to quadratic programming, many iterative and greedy algorithms are available for solving 
the system Eq. fT2TT) . 

Let us first review a basic notion in CS which provides a performance guarantee for BP. 
We say a matrix A £ C mxN satisfies the restricted isometry property (RIP) if 



(1 - 5) \\Zf 2 < \\AZ\\l <(l + 5) \\Z\\l , 5 £ (0, 1) 



(23) 



holds for all s-sparse Z £ C N . The smallest constant satisfying Eq. (T23|) is called the restricted 
isometry constant (RIC) of order s and denoted by S s . 

The following theorem says that the random Fourier matrix satisfies RIP if m is sufficiently 
large. 

Theorem. [25] Let £j £ [0, 1], j — 1, 2, . . . , m be independent uniform random variables. If 



m . „_ 2 ^ 2 . ... 1 



> CtT^s In 2 s In iV In- , rj e (0,1) 
mm 7] 



(24) 



for some universal constant C and sparsity level s, then the restricted isometry constant of 
the random Fourier measurement matrix with 



.4 



e 2 ™^,n 



m 



-N/2,..., N/2 



(25) 



satisfies 5 S < 5 with probability at least 1 — r/. 

Denote X s to be the best s-term approximation of the solution X, and let X be the 
solution of BP Eq.fl22J. 

Theorem. [B] Let A satisfy the RIP with 



5 2s < V2 - 1 



and X be the solution to BP. Then 



X -X 



^Co^WXs-XW. + de , 
\/s 



X -X 



<C r ||X s -X|| 1 + Cie 



(26) 



(27) 



for some constants Co, C\ independent of X. 
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Once the estimate X is obtained from BP, we reconstruct u n by 



u n = -^e- i%2 °X n . (28) 



The problem with Eg. ( 128]) is that the evanescent modes yield exponentially large factor 



n 

cos u + — 

k 



> 1. (29) 



For n sufficiently large, this can magnify the error in X n and produce undesirable result in 
u n . This observation also shows that X may be much more compressible than {u n }. 

A simple remedy would be to apply the hard thresholding by restricting the identity 
Eq.( )28|) up to n sufficiently small and setting the rest of u n zero for \n\ > no- Let us 
now give a rough estimate for the number of modes that should be preserved by the hard 
thresholding rule. 

We define the stably recoverable evanescent modes to be those modes satisfying Eq.( l29p 
and 

k\/3 n \z <C e (30) 
for some constant C e (in [2], C e = 2tt). On the other hand, 



/3 n = ^cos 2 # + 2-cos# + — -1>U-1. (31) 
Hence the stably recoverable modes necessarily satisfy 

l)z <k\P n \z <C e (32) 



n\ 

T 



or equivalently 



C 

\n\ < n = — + k (33) 



which is a rough characterization of the stably recoverable (evanescent) modes. We see that 
no increases as k increases or zq small. 

Summing up the previous analysis we conclude the recoverability of the scattering ampli- 
tude {u n } by the following theorem: 

Theorem. Let z$ > /i max be fixed and let Xj, j = 1,2, ... ,m be i.i.d uniform random 
variables in [— ir, 7r). Let uq = ^ + k for some positive constant C e > 0. Let X, X s be the 
BP solution and the best s-term approximated solution of the system Eq.(|2TT) respectively, 
and assume 

> C— sin 2 sin N In- (34) 



In m 5 2 rj 
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for some universal constant C and n G (0, 1). Let u 
where u n is given by Eq. fj28j) . Then one can reconstruct the solution u with 



V^— no 1 ■ ■ ■ 1 Ujiq ) , u 



no) • • • ) U"n,oj 



\u — u\\ 2 < 



Cn 



(35) 



for some constants Co, C\ with probability at least 1 — 77. 



Proof. Without loss of generality, we prove for the case that Xj, j = 1,2, ... ,m, are i.i.d 
uniform random variables in [0,27r), and consequently the matrix A, defined in Eq. (12 1 [) . is 
the random Fourier measurement Eq. fl25|) where £j are i.i.d uniform random variables in 
[0, 1). It is equivalent to the case where Xj are i.i.d uniform random variables in [— ir, ir): One 
can write Xj = 2ir(£j — |) and the sensing matrix A is then 



>27T(^-|) 



3 27rin§ 



m 



m 



By combining the factor (— l) n into X n and writing W n 



■l) n X n , n 



(36) 



-N/2,..., N/2-1, 



we have 



W-W 



X -X 



\W S - W\\ l = \\X S - X\\ v where W and W s are defined 



in the same manner. 

Under the assumption of the matrix A, we have the estimate 



X -X 



1 



IX.-XUi + Cie 



(37) 



for a desired sparsity level s for some constants Co, C\ with probability at least 1 — rj. On 
the other hand, 



X-X 



n 

E 

ri = — no 



\u„ e ifc/3n20 \/m - u n e ik/3nZ ° v? 



neA 



(38) 



where A = {—N/2, . . . , —no + 1, no + I, N/2 — 1}. Moreover, for |n| < no we have < 



-c, 



e I < | e ifc/3n^o < 1 which ei 



gives 



X-X 



> m J e Ce 1 2 I u n — u n 1 2 + 



"0 



^:II W - U ll2 



(39) 



n=— no 



li- 



no) 



. . . ,w no ) and u = (uZ no , . . . ,u no ). Combining these inequalities, we have 



where u 

the estimate that one can reconstruct u with 



Co^WXs-X^ + de . 



w — u L < 



for some constants C , C\ with probability at least 1 — n. 



(40) 



□ 
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5. Compressibility of the angular spectrum 

Let us now analyze the compressibility of coefficients {u n }. We present a heuristic argument 
suggesting that the angular spectrum of the scattered field is sparse for shallow corrugations. 
Assuming the validity of the Rayleigh Hypothesis we have 

- u inc (x, h(x)) = u(x, h(x)) = u n e ik{a " x+ ^ h{x)) , (41) 

or equivalently 

_ e -ifc/i(x)sm6» _ u n e ik ^ nh ^ e inz (42) 

For sufficiently flat and smooth surface h the nearly normal incidence 9 ~ | tends to produce 
nearly specular diffracted wave [27] and hence {u n } is concentrated at n = 0. This observation 
suggests that it may be reasonable to approximate the outgoing wavevector kf3 n by the 
negative incoming wavevector k(3 , or equivalently, to replace f3 n by f3 = sm9. With this 
approximation, we have 

p — ikh(x) sin(0) 

Y. u ^ nX « ^ oh(x) = - 1 + 2i^(x)/3 + 0(k 2 \h\ 2 ), (43) 

n&Z 

provided that the depth of the corrugation is small compared to the wavelength. Hence, we 
have 



~ A 



" = ° (44) 
2ikh n f3 n^O 



which is sparse by the sparseness assumption on {h n }. Let W = y/m(v n e lkl3 " Zo ) . In view of 
Eq.f lT8|) we have the estimate 

N/2-1 

WX-XsW^WX-WW, = Yl V^\^ nZ0 \ K-v n \ < V^e u (45) 

n=-N/2 

where e u = Yln=-N/2 \ Un ~ v n\- The subsequent numerical simulation shows that u n and v n 
given in Eq. ()44l) are indeed close to each other when the Rayleigh hypothesis is valid. 

6. Numerical simulation 

6. A. Data synthesis 

We compute the scattered field u(x, z ) by the boundary integral method [UJ|2T]. The scat- 
tered wave can be represented by the Brakhage- Werner type ansatz, i.e. the representation 
via mixed single-layer (S) and double-layer (K) potentials 

u = (K - ir)S)<p on fi r (46) 



with a mixed layer density ip for a constant 77 > which can be adjusted to improve the 
condition number of the system. Explicitly, we can write 

u(x,z) (47) 
-^®((x,z),(x,h(x'))) -ir ] <f>((x,z),(x,h(x'))) S J ijj{x') ■ yjl + h 2 (x') da/. 

Taking the limit z — > h(x) + and using the properties of single and double layer potentials, 
we obtain the boundary integral equation [9|[TTj 



u™(x,h(x)) = U{x)+ r (J^$((x,h(x)),(x',h(x'))) 



ir]$({x,h{x)),(x',h(x'))) ) ip{x')\Jl + h 2 {x!)&x! (48) 



(see Appendix). Note that the integral in Eq. (l48[) has weakly singular kernel, and the integral 
exists as an improper integral since the periodic Green's function 

$(r,r') = l -^e 2nmkcose H^\k |r - r 1 - 2tto(1, 0)|), x, x' E [— tt, tt) (49) 

raSZ 

has the same singularity as Hq(x) ~ 1 + — (ln| + 7) where 7 0.5772 is the Euler- 
Mascheroni constant. Moreover, 

±H?\k \r - r'\) = -kHf\k \r - r>\) ^±^ (50) 

converges to a finite limit (a curvature-like term w.r.t. the boundary) as r — > r' [10] implying 



With i/j solved from Eq. fl48p and the Sommerfeld integral representation 



the boundedness of on T. 



H^(k\r\) = ij e ifc d^+-)^, (51) 



where 



p = \y^i \«\<] (52) 

[ Wo 2 - 1, \a\ > 1 
we obtain from Eq. (j4"7"|) the outgoing wave expansion for the scattered field 

u(x, z) = J2e ikianX+M ( ^ Q-^nx'+Mx')) 3 n ( x ')^( x ') dx >\ j z > /, max (53) 

with 



g n (x') = k - kh(x')^ + T^-V 1 + h 2 (x'). (54) 

Pn Pn 

9 



Comparing Eq. (1531) with Eq. (|T2|) we arrive at the expression 

u n =^j: n e-< a - x ' + ^ hix ''>)g n (x')i;(x')dx' (55) 

relating the angular spectrum of the scattered field to the mixed layer density ip. 

Eq. (l48l) and Eq. (!55|) motivates the following iterative reconstruction scheme. Given h^ m \ 
m = 1, 2, 3, . . . first solve for if)^ from 

-u inc (x,h {m) (x)) = ^ m \x) + J* (J^<S>((x,h( m \x)),(x',h( m \x'))) 

-ir]$((x, /i (m) (x)), (x f , h {m \x')))^ {m \x')^l + \h( m )\ 2 (x') dx(56) 

and then solve for /?,( m+1 ) from 

u n = -i f\M a ^'+^ him) ^)g^+ 1 \ x ')^( x ')dx\ neZ (57) 

47] " «/ — 7T 

g^ +l \x') = k - kti m+1) (x')^ + + \h( m )\ 2 (x'). (58) 

Pn Pn 

Note that both Eq. (1561) and Eq. (1571) are linear equations. 

A natural candidate for the initial guess of the above iteration is the one obtained under 
the Rayleigh hypothesis that the validity of Eq.f l53|) is extended to the region z > h(x). 
Specifically, we extend Eq.( l53~]) all the way to boundary and study the nonlinear equation 
Eq. pTl) . Indeed, this alone produces excellent results for shallow corrugations and will be the 
focus of the following numerical experiments. Scattering and imaging of shallow corrugations 
can also be treated by assuming the Born approximation |12j . 

In our numerical simulations, we set 128 nodes to solve the boundary integral equation 
Eq. (T48!) by the Nystrom method, with 77 = 1. Figure [2] shows two examples of the computed 
scattered field. We define R a (h) = max n 2\nh n \ for a rough metric of the validity of the 
Rayleigh hypothesis. 

6.B. Surface reconstruction 

Solve for h(x) from Eq. (HT|) . we consider the following three algorithms: the first two are 
pointwise matching schemes and the third is a global fitting scheme. 



1. Point-wise, fixed-point iteration for /i(x),Vx G [—it, it). A fixed point iteration algo- 
rithm was introduced in [Tt)J[27] and is described below. The initial condition (x) is 
chosen in the following way. For 9 rs | the angular spectrum {u n } is concentrated at 
n = 0. Substituting (3 n by 1 in Eq.f HT]) yields 

e ikh(x) _ u e ika n x _ _ Q ik (a; cos 9~h(x) sin 6») 

n&Z 
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Fig. 2. The real part of the scattered field induced by the profile h(x) = 
0.2454 sin(x) (left) and the profile h(x) = 0.01 sin(12x) + 0.007 cos(7x) (right) 
with L = 2tt and vertical incident wave. 



One solves Eq. (1591) by the iterative scheme 



x) 



-2ik 

111 (- EneZ U n? 



ik(a n x+(l3 n -l)hW(x)) 



-2ik 



for n — 1, 2, . . . and all x G [— 7r, it). 
2. Newton's method. From Eq.( l4T|) . for each x G [—ir,ir) we set 



and solve e(/i, x) = for /i(x) by Newton's method 

e{h®;x) 



h H+i] = h [H 



±e{h^x) 



with initial value Eq.( )60|) . 
3. Nonlinear least squares fitting. Let 



F(a) 



^ \e(^2a n (f) n (xj),x 



2 



(60) 
(61) 



(62) 



(63) 



(64) 



where e(h; Xj) is defined in Eq.( )62l) and a — (ai, 02, . . . ) is the vector of the coefficients 
of an expansion of h corresponding to a frame {4> n }, i-e., h = ^2 n a n (j) n . Consider 
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minimizing the nonlinear least square 

minF(a). (65) 

a 

The basis function are chosen to be sin(nx) and cos(nx) for n £ II C N, where the 
index set II contains those indices n such that |iT n | are relatively large. The Matlab 
subroutine lsqnonlin is applied, which is based on the subspace trust region method. 

6. C. Examples 

In the following examples we apply vertical incident wave, 6 = |, with the wave number 
k = 3.2 (i.e., the wavelength A ~ 1.9635). After synthesizing u(xj, z ) for j = 1,2, ... ,m, 1% 
additive noise, with respect to ||it(-, Zo)\\2/\^n, is added to the measured scattered field u. 
The profile functions h are L = 2tt periodic, defined in one period [— ir,ir), and periodically 
extended into M. The bound for the exponential factor is set to be C e = log(25). Yalll 
[29] algorithm, a Basis Pursuit solver, is applied to solve Eq. fl20|) for vector X. To avoid 
exponential amplification of small components of X, we apply a threshold level r = 10% • 
max n7 io \X n \ and filter out the components below r, then compute u n from Eq. (l2"H|) . 

In Figures [MHl the right panels show the exact profile h(x) (black solid line), and the 
reconstruction under the three algorithms: Newton's method "Newton", fixed point iteration 
"Fixed pt iter", and nonlinear least squares fitting "NLS fit". The length of the black 
strip on the top of each plot indicates the wavelength, and its height indicates the vertical 
coordinate Zq of the sampling points. The left and middle panels show the real (left) and 
imaginary parts (middle) the angular spectrum u n (blue crosses), the estimated angular 
spectrum u n (red dots), and the theoretical estimate v n (green circles). Note the different 
scales for the real and imaginary parts in Figures [HE where R a is relatively small. This is no 
longer the case in Figures [7] and |8] for which the Rayleigh hypothesis is known to be false. 

Figures [3] and H] show the results for profiles h(x) with sparse Fourier coefficients. The 
prediction v n captures well the dominant component dt[uo] and so does the sparse recon- 
struction u n the other significant components of the angular spectrum. For reconstruction 
(right panels) the nonlinear least squares is the best performer while the pointwise iterative 
methods may produce visible undershoots at the peaks and troughs. 

For the Gaussian profile (FigEJ) and subwavelength double-peaks (Fig|6]), again v n captures 
well the dominant component and so does the sparse reconstruction u n most other 

significant components of the angular spectrum. The angular spectrum for the latter case 
occupies a wider range of modes than the former case since the two peaks are sharper than 
the Gaussian. As a consequence, the reconstruction is more accurate in the former case. For 
the latter case, all three reconstructions undershoot the peaks and produce fluctuations at 
the flat part of the profile. 

Figures [7] and E] are the results for simple sinusoids when the Rayleigh hypothesis is known 
to fail (ab > 0.448). The failure of the Rayleigh hypothesis manifests in the broadening of the 
support of the angular spectrum. Furthermore, the imaginary part of the angular spectrum 
is order of magnitudes larger than those in Figures [3J16J As a result, the angular spectrum 
{u n } is less compressible and not well recovered by the compressed sensing techniques. In 
both cases, the simple prediction v n fails to capture the dominant components of the angular 
spectrum. 



12 



Angular Spectrum 



Reconstruction, k=3.2, e.=0.5ji, z Q =0.055 
Ra=0.13, e =0.0479736, m=64, n =62 




Fig. 3. The profile of 5 Fourier modes h(x) = 0.0l(^ sin((l + 3p)x)) and 
the reconstructions. 



Angular Spectrum 



Reconstruction, k=3.2, 0=0.5?!, z Q =0.02 
Ra=0.12, e =0.00688774, m=64, n„=165 



-100 



-100 




Fig. 4. Profile of two Fourier modes h(x) = 0.01 sin(12x) + 0.007 cos(7x) and 
the reconstructions. 
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Angular Spectrum 



Reconstruction, k=3.2, EMJ.&E, z =0.02 
Ra=0.00482198, e =0.00903647, m=64, n n =165 




***** 



- Exact 
Newton 
Fixed pt iter 
NLSfit 



Fig. 5. Periodized Gaussian h(x) = b(e (ax)2 - E ^ L ) ■ X[-o.97t,o.9vt] (%) , a = 2, 
b = 0.01 and the reconstructions. Here x is a smoothed indicator function. 



0.05 

■ 



Angular Spectrum 



-200 -100 100 200 



Reconstruction. k=3.2, 0=0.5:1, z Q =0.012 
Ra=0.01 11735, e =0.0277698, m=64, n n =272 



-200 -100 100 200 



-200 -100 100 200 




Fig. 6. Double subwavelength peaks h(x) = &(C(a(x— |))+C(a(x+|))), a = 2.5, 
b = 0.01, = exp (l — T^rry) X(-i,i)( x ) + c o and the reconstructions. Here 
the constant c is chosen such that Co = 0. 
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Reconstruction, k=3.2, 6==0.£m, z q =0.736311 




Fig. 7. h(x) = 0.491 cos(x) 



Reconstruction, k=3.2, ©.=0.571, z Q =0.883573 
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Fig. 8. h(x) = 0.589 cos(x). 
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Nevertheless, the nonlinear least squares fitting provides an accurate reconstruction of the 
profile in both cases. The Newton iteration converges in Figure [7] but fails near the peaks 
and troughs in Figure M while the fixed point iteration fails to converges near the peaks and 
troughs in both cases. When ah is further increased (to, e.g. 0.736), then all three methods 
fail to recover the profile. 

7. Conclusion 

We have proposed a compressed sensing scheme for near-field imaging of corrugations of 
relative sparse Fourier components. The scheme employs random sparse measurement of 
near field to recover the angular spectrum of the scattered field. We have shown heuristically 
and numerically that under the Rayleigh hypothesis the angular spectrum is indeed sparse 
or compressible and amenable to compressed sensing techniques. 

We then develop iteration schemes for recovering the surface profile from the angular 
spectrum. Specifically, under the Rayleigh hypothesis we have tested three iterative schemes. 
The nonlinear least squares in the Fourier basis has the best performance among the three 
and produces accurate reconstructions even when the Rayleigh hypothesis is known to be 
false. 

The full iteration scheme Eq. fl56|) -Eq.f l57|) beyond the limitation of the Rayleigh hypothesis 
will require non-sparse measurements for the angular spectrum data and will be studied 
elsewhere. 
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A. Derivation of the boundary integral equation Eq. (148ft 

The term \ip{x) in Eq.( l48p arises due to the jump discontinuity for the double layer potential 
across the boundary, whereas the single layer potential is continuous. More specifically, let 

u s (r) = J <$>(r,r')ij;s(r')dS(r') (66) 

—$(r,r')i> D (r')dS(r') (67) 

be the single and double layer potentials respectively for r = (x, z) G M 2 \ T. Furthermore 
we denote r ± = Vq ± pz/(ro) for some small p > and Vq e T (assuming that the boundary 
is of class C 2 so the representation of r ± is unique for r near the boundary). Clearly 



\imus(r + ) = \imus(r ) = us(r ). (68) 

p-s>0 p->0 



On the other hand, write 



u D (r + ) = Mro) J ^M r+ ' r') dS(r') + v(r + ) (69) 
—<S> ( r -,r')dS(r')+v(r-) (70) 



so that 



J ^(r* r') {Mr') - Mro)) dS(r') 

+ Mro) jT (^7 $ (^> O - ^7*0^, r')\ dS(r') 



(71) 



where <J>o is the Green's function for Laplace equation. It is easy to see that integral Eq. (l71j) 
is continuous in the neighborhood of p = 0. 
The jump condition 

lim (u D (r+) - u D (r~)) = M^) (72) 
now follows from the calculation 



/ tJUo^, r') dS(r') = \f -^-M r± > O dS ( r ' 
Jr ov' 2 J dBp {r ) ov> 



(73) 



/ ±ldS(rO— >±J, P^O (74) 



by applying the divergence theorem, integrating over the circle B p (r ) of radius p and shrink- 
ing radius p to 0. 
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